library(ggplot2)
library(reshape)
library(rjags)

data<-read.csv("all-votes.csv",header=T,as.is=TRUE)
## The data must be ordered with cases lacking directionality to the top
my.order<-order(data$direct1,na.last=FALSE)
data<-data[my.order,]

## Now we must (get ready to) tell jags which cases lack directionality
case<-as.numeric(factor(data$casenum,levels=unique(data$casenum)))
missing.direction<-max(case[is.na(data$direct1)])
has.direction<-missing.direction+1


## And also about the areas and the directions
area<-as.numeric(factor(data$area))
direction<-tapply(as.numeric(data$direct1==1),case,unique)
area<-tapply(area,case,unique)
no.areas<-max(area)
## Prepare the list for Jags
forJags<-list(y=data$vote,
	case=case,
	no.cases=max(case),
	missing.direction=missing.direction,
	has.direction=has.direction,
	area=area,
	no.areas=no.areas,
	direction=direction,
	justice=as.numeric(factor(data$judge)),
	nObs=nrow(data))

## inits

data.c<-cast(data,casenum~judge,value="vote")
judge.names<-names(data.c)[2:ncol(data.c)]
data.c<-as.matrix(data.c[,2:ncol(data.c)])
data.c[is.na(data.c)]<-mean(data.c,na.rm=T)
judge.inits<-as.vector(pscl:::x.startvalues(t(data.c),d=1))
names(judge.inits)<-judge.names

case.inits<-pscl:::b.startvalues(t(data.c),judge.inits,d=1)

set.seed(1982)

system.time(foo<-jags.model("hier.jags",
    n.adapt=1e3, ## 
    data=forJags,
	inits=list(theta=judge.inits,
	beta=case.inits[,1],
	alpha=case.inits[,2])
)
)

	max.iter<-1.25e6
	max.burnin<-.75e6

my.thin<-max.iter/4000

system.time(update(foo,n.iter=max.burnin)) ## 5984s
out<-coda.samples(model=foo, 
    variable.names=c("theta","beta","alpha","lm.coef","lm.intercept"),
    n.iter=max.iter,thin=my.thin) 

source("diagnostics.r",echo=TRUE)

	save.image("hier.RData")

